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Abstract 

We introduce the new concept of an EBV to assess the sensitivity of 
model outputs to changes in initial conditions for weather forecasting. 
The new algorithm, which we call the Ensemble Bred Vector or EBV, 
is based on collective dynamics in essential ways. As such, it keeps 
important geometric features which are lost in the earlier bred vector 
algorithm (B V) . By construction, the EBV algorithm produces one or 
more dominant vectors, and is less prone to spurious results than the 
BV algorithm. It retains the attractive features of the BV with regard 
to being able to handle legacy codes, with minimal additional coding. 

We investigate the performance of EBV, comparing it to the BV 
algorithm as well as the finite-time Lyapunov Vectors. With the help 
of a continuous-time adaptation of these algorithms, we give a theo- 
retical justification to the observed fact that the vectors produced by 
BV, EBV, and the finite-time Lyapunov vectors are similar for small 
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amplitudes. The continuum theory is estabhshes the relationship be- 
tween the two algorithms and general directional derivatives. 

Numerical comparisons of BV and EBV for the 3-equation Lorenz 
model and for a forced, dissipative partial differential equation of 
Cahn-Hilliard type that arises in modeling the thermohaline circu- 
lation, demonstrate that the EBV yields a size-ordered description of 
the perturbation field, and is more robust than the BV in the higher 
nonlinear regime. The EBV yields insight into the fractal structure 
of the Lorenz attractor, and of the inertial manifold for the Cahn- 
Hilliard-type partial differential equation. 

Keywords: Bred vectors, Lyapunov vectors, sensitivity, dynamic stability, 
Cahn-Hilliard, Lorenz. 
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1 Introduction 

Central to weather prediction is the analysis of the sensitivity of a physical or 
computer-coded model to initial conditions. Model sensitivity to parameters 
is also important in model inter-comparison. One studies such sensitivity in 
order to obtain a better understanding of the role played by these parameters 
in model outcomes. 

Sensitivity and predictability are often intertwined in the context of weather 



prediction and have been the subject of extensive research (see Buizza et al. 



(1993) and references contained therein.) These are not exclusively weather- 
related issues and thus geophysical fluid dynamics will often mine other phys- 
ical, computational and mathematical disciplines, for ideas with which to 
assess dynamic sensitivity. Practical sensitivity methodologies must con- 
tend with the evolution and dynamics of highly coupled, complex, high- 
dimensional systems, riddled with subscale parameterizations and empirical 
relations, which are the norm in large-scale climate and meteorology models. 

A tool used in the study of sensitivity analysis is the Bred Vector (BV) 
algorithm. It is proposed for use in forward sensitivity of weather and climate 
models. While in Subsection 2.1, we present a brief survey of some of the 
applications of this algorithm in various sensitivity analyses, in this article, 
we will focus on the issue of the maximal growth of errors due to small 
changes in the initial conditions. 

The concept of the BV algorithm we use is based on the theory first 
introduced in Toth and Kalnay (1993). In addition to the BV notion, we 
present here a new variant, which we call the Ensemble Bred Vector (EBV) 
algorithm. The definitions of both the BV and EBV algorithms are presented 
in Section 121 



In the BV algorithm, one follows an initial condition of the time-discrete 
nonlinear system, along with cloud, which describes a family of nearby solu- 
tions. (Since this algorithm is used to sample the error space, an ensemble of 
initial perturbations is bred simultaneously.) The perturbations at the initial 
time are fixed with a common small amplitude e. After each cycle, the out- 
come of the perturbations is rescaled to the same amplitude e. For the BV 
algorithm, the rescaling of each perturbation is independent of the others, 
and there is no mechanism to use the rescaling to compare the dynamics of 
nearby perturbations. 

The new variation that we propose here, the EBV algorithm, differs from 
the BV algorithm, in the rescaling rule. In particular, for the EBV those 
perturbations that are not the same size as the largest perturbation, play a 
reduced role after the rescaling. Thus the rescaling used in the EBV algo- 
rithm serves us better in separating various levels of the dominant dynamics. 
In short, the EBV algorithm offers better insight into the relative behavior 
of nearby trajectories. Therefore, even when the initial perturbations of the 
two algorithms are based on the same cloud, the EBV algorithm is linked to 
the ensemble dynamics of the underlying non-linear model, hence its name. 
We will make use of both the BV and the EBV algorithms. In fact, one of 
our major goals is to present an in-depth comparison of the two algorithms, 
as they are used, or can be used, in describing the underlying dynamics of 
the model. 

We will show that in some metrics, the BV and EBV algorithms are 
comparable, with the EBV being more accurate and faster, see Table [T] and 
Figures |6| [7| [8} for example. For other issues, especially those involving 
the longtime dynamics within the global attractor, the BV algorithm has a 
shortcoming, which limits its use (see Section 4). The EBV, instead, leads 
to useful and interesting insight into the dynamics of the model, as is shown 
in the three Figures |4| [T], and |2} 

This brings up the question: For a given model of sensitivity with respect 
to initial conditions, how does one determine the direction vector that results 
in the maximal increase in the error due to a small perturbation in the given 
direction initially? This is where the Lyapunov vectors enter the scene. What 
one needs is a red vector, which is the Lyapunov vector xq with ||xo|| = 1 
and with the property that the corresponding strong Lyapunov exponent Ai 
is the maximal Lyapunov exponent for the model. (See Subsection 2.5 for the 
definition and more details. One should note that the Lyapunov exponents 
for the model require integration over < t < oo, or over the real line M.) 

We will use either the EBV or the BV algorithm to approximate the 
solutions of the tangent linear equation. In Section 3 we show that either 
algorithm is a good approximation. For these algorithms one can integrate 
only over a finite interval < t < T, where < T < oo. However, it 
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is only under exceptional circumstances, e.g., time-periodic or autonomous 
problems, that a finite-time integration will approximate well an infinite time 
average. For any hope for success in using a finite-time integration, we require 
that the model problem satisfy two properties: 

(1) There is an attractor 21 for the solutions, 

(2) The initial condition is chosen to be on, or near, 21. 

What can one expect with such a finite-time approximation, when the 
initial condtions are near the attractor 21? In terms of the calculated time, 
one expects first to be in a transient state. Then after a while, one hopes to 
get some meaningful information about the longtime dynamics of the model. 
We include in this manuscript several studies of such approximations. 

It is very important to note that, in order to better understand the maxi- 
mal growth of errors due to small changes in the initial conditions, one needs 
to exploit the dynamical information contained in the attractor of the model. 
In particular, one needs to complete two steps: 

• Step 1: One needs to locate the red vector xq and the associate Lya- 
punov exponent Ai. This then determines the red vector solution 
R{xo,t), for — oo < t < oo, for the model equation. 

• Step 2: One must find a good approximation of the time evolution 
R{xQ,t) on an appropriate finite-time interval, < t < T < oo. 



Once the red vector xq is known, then e xq is the initial condition for the 
bred vector sequence, see Subsection 2.2. Due to the results derived here in 
Section 3, either EBV{t), or BV{t), is a good finite-time approximation of 
R{xQ,t), for < t < T < oo. This takes care of Step 2. Consequently, the 
problem boils down to a search for the red vector, which is addressed below 
in Subsection 2^ Fortunately for us, there is wealth of related mathematical 
information in the 1987 manuscript of R. A. Johnson, K. J. Palmer, and G. 
R. Sell, see Johnson et al. (1987). (We will refer to this paper as the "JPS87" 
in the sequel.) As we shall show, by using the JPS87, we are able to describe 
the mathematical process of finding the red vector. By using this citation, 
with the theory of the EBV{t) algorithm, this leads to a good solution for 
our sensitivity problem. 



It was observed by Toth & Kalnay op. cit. (see also Toth and Kalnay 



(1997)) in several experimental runs that BVs resemble the leading finite- 
time Lyapunov vectors. In order to make a more quantitative comparison 
between Lyapunov vectors and BVs, in Section |3] we study the Continuum 
Limits (as the basic time-step size for rescaling goes to 0) of the BV and 
the EBV algorithms and show direct connections between these limits and 
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specific solutions of tlie continuous-time tangent linear equations. Section 
[3] also contains a further discussion of some of the desirable and interesting 
features of the new EBV. For instance, there is a natural ordering of the 
ensemble members of an EBV, and as we will show, it is possible to observe 
perturbations with smaller sizes than the dominant one, but with very strong 
growth, see the spear-like behavior in Figure |4j 

We consider two models to exemplify the features of the new EBV. The 
first is the familiar Lorenz63 model introduced by Lorenz (1963). It has a 
well-known global attractor. The second is a nonlinear forced and dissipative 
partial differential equation of the Cahn-Hilliard type. This equation is a 
variation of a model proposed by Cessi and Young (1992) of the oceanic 
thermohaline circulation. We will denote the equation associated with this 
model as the CY92 in this study. (In fact, we impose periodic boundary 
conditions instead of the more physical zero-flux, zero-stress conditions at 
the poles.) 

It turns out, it is a good example of the typical climate- related model 
dynamics. However, the CY92 is special, since it has an inertial manifold. 
Consequently, the longtime dynamics of this partial differential equation is 
completely contained in the attractor of a finite dimensional ordinary differ- 
ential equation. In applications, the BV algorithm and its variants are in 
fact discrete-time algorithms based on finite-dimensional approximations of 
weather models, obtained either by mode projection as in the Lorenz (1963) 
system, or by spatial discretization of partial differential equations (PDEs), 
as in the Cessi- Young (CY92) model, which will be described below. As we 
will show the EBV yields insights into the structure of the attractor of the 
Lorenz63, and of the inertial manifold to the CY92. 

Our numerical examples also highlight that the BV algorithm is sensitive 
to the amplitude and frequency content of the initial perturbation. In con- 
trast, the outcome of the EBV algorithm shows a clear hierarchy among its 
members, and the first few members already generate an unambiguous char- 
acterization of the perturbation field at both large and small amplitudes. In 
the nonlinear regime, however, the EBV will be shown in Section |4]to be less 
likely to produce spurious results than the BV. In Section |5] we will address 
implementation issues of the EBV. 



2 BV and EBV Algorithms; Finite-Time Lya- 
punov Vectors 

In this section we present the definitions and methodology for computing 
the two Bred Vector algorithms, the BV and the EBV. We also review the 
basic theory of Lyapunov Vectors, Lyapunov exponents, and their finite-time 
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counterparts. We will compare these different tools below. However, before 
doing this, we include here a brief survey of some of the applications and 
theoretical issues that have been noted in the implementation of the BV 
algorithm. 



2.1 Brief Survey of Applications of the BV Algorithm 



For the BV algorithm, we use the one originally proposed by Toth and Kalnay 



(1993). BV is purely algorithmic. It is "equation- free" and thus with addi- 



tional minimal computer coding it can handle legacy code representing even 
extremely complex models. Most alternatives for obtaining estimates of for- 
ward sensitivity will involve non-trivial additional coding. For example, in 
order to obtain the finite-time Lyapunov vectors and exponents, one needs 
to derive and make use of the tangent linear model. (From the beginning, 
it was realized that there existed some close connections between the BV 
algorithm and the tangent linear equations. As we will show in Section 3, 
there is a rigorous mathematical foundation for these connections.) Singular 
value decomposition methods, which can offer complementary information 
to Lyapunov-vector inspired methods, also require a tangent linear model. 



Deremble et al. (2009) use this approach to study regime predictability in 



some reduced weather models. Other examples are: Buizza et al. (1993) and 



Palmer et al. (1998). In another direction, Wolfe and Samelson (2007) pro- 



pose the use of the MET and the finite-time singular vectors to approximate 
the Lyapunov vectors. 

The BV algorithm is a finite-time, forward sensitivity methodology which, 
in addition to being useful in characterizing model sensitivity to initial condi- 
tions, has been proposed as a means to produce a reduced-rank representation 
of the background error in data assimilation and forecast error- covariance ap- 
proximations 



see 



Corazza et al. (2003), for example 



Several articles in the literature have addressed applications of the BV 



algorithm in weather modeling. See, in particular Toth and Kalnay (1993) 



Toth and Kalnay (1997), Kalnay (2003), and reference therein. For a compar- 



ison of the BV algorithm and other methods, such as Monte-Carlo perturbed 



observations, we refer for example to Cheung (2001); Gneiting and Raftery 



(2005); Hansen and Smith (2000); Wei and Toth (2003). For an applica 



tion of the BV algorithm to ensemble Kalman filters see Wang and Bishop 



(2003). Primo et al. (2008) and Hallerberg et al. (2010) have treated applica- 



tions based on variations in the algorithm, where for example, the rescaling 
is done by using a geometric mean. 

As already discussed, several alternatives to the BV algorithm have been 
proposed and employed in the literature. BVs have been viewed as non-linear 
analogs of finite-time Lyapunov vectors. Similarly nonlinear analogs of sin- 
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gular vectors have been proposed, for instance conditional nonlinear optimal 



perturbations proposed by Mu and Jiang (2008) and non- linear singular vec- 



tors proposed by Riviere et al. (2008 ), although they entail a computationally 
expensive optimization. 

Both the BV and EBV algorithms arise in the time-discretization of a 
continuous-time dynamical system. We consider the following initial value 
problem 

I = (1) 

z/(0) = yo, 

where t represent time and G = G{y) is a map that has at least a bounded 
gradient. Since our main applications involve autonomous differential equa- 
tions, we assume that G does not explicitly depend on time. The basic theory 
we present here has a routine extension to non-autonomous problems. 

The solution vector y = y{t) can live in a finite- or infinite-dimensional 
normed linear space. In the former case, ([T]) is an (autonomous) system of or- 
dinary differential equations, while in the latter case, ([T]) is an (autonomous) 
system of partial differential equations, modeling a time dependent, spatially 
extended system. For systems of partial differential equations, we assume 
that either periodic boundary conditions or non-flux boundary conditions 
are prescribed. See for example. Sell and You (2002). The use of other 
boundary conditions may lead to a related theory, but we do not address the 
issue here. 

If the system contains evolution partial differential equations, the system, 
along with the boundary conditions, are discretized in space or projected onto 
a finite-dimensional space compatible with the boundary conditions. Conse- 
quently, we usually assume that ([T]) is a system of ordinary differential equa- 
tions of dimension K, which may be large. Since most large-scale weather 
and climate circulation models presently use explicit-in-time integrators, we 
will focus on numerical models of this type. 



2.2 Bred Vector Algorithms: BV and EBV 

We let y = y{t) denote a given (continuous time) solution of ([!]). We then 
turn to an approximate solution Yn = Y{tn), which is defined on the time 
grid: t„, for n = 0, 1, 2, ■ ■ ■ , where = tn + Stn, for = 0, 1, 2, . . .. We set 
to = 0. We assume that 5t„ = 6t is positive and small, and that it does not 
depend on n, for n > 0. 

For the autonomous case, the initial value problem, which is approximated 
using an explicit numerical integration scheme, leads to consideration of the 
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difference equation 



n+l — 



M(r„,5t), n = 0,l,2,.., 



(2) 



Yo = Y{0). 



where Y(tn) is a solution of the discrete problem ([2]), and it may be viewed as 
an approximation of the continuous-time solution y{t), at t = t„. Likewise, 
the initial condition Yq is an approximation of the initial condition yo. The 
points Yn are in the i^-dimensional Euclidean space R-^ and M = M{Y,6t) 
is the discrete-time solution operator on generated by the ordinary dif- 
ferential equation ([!]). As noted above, we assume that the discrete-time 
problem ^ has an attractor 21, and that Yq is on, or near 21. 

At this point it is convenient to introduce the related concepts of a Cloud 
(at t = 0) and a Family of Initial Perturbations. A Cloud (at t = 0) is a 
family of tangent vectors Syo{L) to Yq in that depends on l, where l El 
and I is a finite index set. The main requirement we impose is that 



where e is small, positive, and fixed. The collection of all terms (Yq, Syo^i)) in 
{Yq} X M^, for i G I, is called a Family of Initial Perturbations. Notice that 
this collection lies on a sphere of radius e in {Yq} x M.^ with center at [Yq, 0). 
Both the Cloud and the perturbations evolve in time, via the BV or the EBV 
algorithms, which are defined below. As we now note, the definitions of the 
two algorithms differ only in the rescaling rule. 

We begin by recalling the BV algorithm as given by Toth&Kalnay. For 
n > 0, we assume that the base point Yn and the perturbation vector 5F„(i) 
are known. For the (n + 1)*^* step we use: 

1. Yn+i denotes the {n + base point, and it is determined by ([2]); 

2. 6yn+i{L), the (n + 1)*^* perturbation vector, is given by 



where Rn+i is a rescaling rule. 

The time evolution of the Cloud is BV(t„) = 53^„(t), where l El and n > 0. 

One rescaling rule, proposed by Toth and Kalnay, op. cit., consists of 
rescaling the n^^ perturbation vector to the previous one by 



\\Syn+i{i)\\ = Rn+i\\SYn+i{L)\\ = \\Syo{t)\\ = e, foiLEl,n> 0. (5) 



11^^0(011 



= e 



for all L eI, 



(5r„+i(0 = M{Yn + 6y^{L),6t)-M{Y^,6t) 

53^n+l(0 = Rn+1 6Yn+l{L), 



(3) 
(4) 



Equivalently, 



Rn+l 



5K„+i(0|| 



ii^yo(oii 



(6) 
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so that Rn+i = Rn+i{Yo) depends on the initial base point Yq, as well as 
the initial perturbation vector 6yQ{L). An alternate rescaling rule consists in 
rescaling periodically, at t^k, where k is an integer k > 2, and m = 1, 2, ■ ■ ■ . 
In this case, one uses the rule ^ when n = mk, and 

Rn+iiYo) = 1, for mk <n < {m + l)k. (7) 

We do not use this alternate rule in this paper. (For a discussion of rescaling 



time and regime predictability in some reduced models, see Deremble et al 



(2009) and references therein). 

For the Ensemble Bred Vector algorithm, instead of using Rn+i as in 
equation ([s]), we use a uniform scaling -R™"i, which is the same for all t G I. 
In particular, we replace equation ^ with 

5F„+i(0 = M(F„ + 53^„(0, 6t) - M{Y^, 6t) 



for all i G I, where 



"D mm 



max(\\6Yn+iU) 



(9) 



Similarly, when ^ and ^ hold, we use EBV(t„) = Synii), for i G I and 
n > 0, to denote the time evolution of the Cloud. (Alternatively, a periodic 
rescaling rule utilizing ([t]) above can be employed.) The time step used to 
compute the base trajectory and the time intervals between normalizations 
need not be the same. This is the case for both the BV algorithm as well as 
EBV. 

The crucial difference between the BV and EBV algorithms is that, even 
when the BV is run concurrently over an ensemble of initial data, the out- 
come of the algorithm for each given datum does not depend on the other 
members of the ensemble. In contrast, the evolution of the ensemble mem- 
bers is interdependent in the EBV. Nevertheless, by construction, the EBV 
should exhibit similar behavior to the BV at small amplitudes. Indeed, one 
of the design principles for the EBV algorithm is to reinforce this aspect by 
providing it with a built-in acceleration mechanism. 

Both the BV algorithm and EBV outcomes, on the other hand depend 
on the choice of vector norm used to define the rescaling rule in either ^ or 
(|9|. While all norms are equivalent in the (finite) i^-dimensional space where 
we seek solutions, in practice the constants appearing in the equivalence 
between different finite-dimensional norms generally strongly depend on the 
dimension K and eventually blow up as k becomes infinite. Hence, the choice 
of norm used can have an impact on the implementability and performance of 



these algorithms. In Riviere et al. (2008), it was suggested that the outcome 
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of the BV algorithm may depend strongly on the choice of norm. This is 
actually a consequence of the non-negligible nonlinear effects in the system. 
In Section |5| we will elaborate further on the issue of norm dependence. 

In addition to the "renormalization time-step" 5t, there is another time 
step, the "integration time-step", which we will denote by At. For example, 
one encounters the new time step when moving from the continuous-time 
problem equation (1) to the discrete-time problem equation (2). For the most 
part, we will treat At and 5t as being equal in the calculations described in 
this article. However, we always require that At < 5t. 



2.3 Lyapunov Vectors 

In Section |4| we will be applying the B V and EB V algorithms to two models 
consisting of systems of (nonlinear) autonomous ODEs of the form ([T]). ( The 
second model arises from the discretization of a PDF.) In each model, there 
is a compact, global attr actor 21, which is a subset of M^, and is invariant for 



the time evolution of the system. (We refer the reader to Chapter 2 in Sell 
and You (2002 ) for more information on attractors and global attractors.) We 
will let 9 denote a typical point in the attractor 21, and we will let ^-t = y{t) 
denote the unique solution of ([T| that satisfies ?/(0) = 6. Since 21 is invariant, 
one has 6^ ■ t e 21, for all t G M. 

In order to study the sensitivity with respect to initial conditions on the 
attractor, the Tangent Linear Model is used, which is defined as 



dtx = A{e ■ t) X, 



(10a) 



x{0) = xo G M^, (10b) 

where A{y) = DG{y) is the Jacobian matrix of G. Hence A = A{6 ■ t) is 
the linearization of ([T| along the solution 9 ■ t. We observe that, even if G 
does not explicitly depend on time, (10a) is generally non-autonomous, since 
6 ■ t changes with time. 

We let U{9,t), denote the solution operator of ( |10a ), which takes the 
initial data to the solution at time t, so that U{6,t)xQ is the solution of 
the initial value problem for (10). Such an operator is well defined by the 
uniqueness of solutions to the problem (10). Uniqueness of solutions also 
readily implies the cocycle identity: 



U{9,T + t) = U{9 ■r,t)U{e,T), 



for all ^ G 21 and all t, t E 



'111 



Next we consider a family of mappings H = H(t), which are defined for 
t G M by the relation 

H(t)(^,Xo) := (^■t,?7(^,t)a;o), for (^, xq) G 21 x M^. (12) 
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We note that 11 (t) maps 21 x R-^ into itself, for each t G 
continuous in {t,9,XQ); it satisfies n(0)(6',a;o) = {9,xo), (i.e. 
identity operator; as well as the evolution property: 



n(r + t) = n(t)n(r) 



for all T, t E 



it is jointly 
n(0) = /, the 

(13) 



By using the discrete-time dynamics, where t and r are restricted to 
satisfy t = n ■ 6t and t = m ■ 6t, the notation and the theory of dynamical 
systems extends readily to the discrete-time problems of interest herein. Note 
that the ^-component of 11 does not depend on the xo-component. Thus 11 
is called a skew product flow. Since 11 is linear in xq, it is sometimes called a 
linear skew product flow. In summary, the Tangent Linear Equation over the 
attractor 21 generates a linear skew product flow. (For more information on 
the theory of skew product flows in the context of non-autonomous dynamics, 



see the multiple works of Backer and Sell, for example: Sacker and Sell (1977) 



Sacker and Sell (1978) and Sacker and Sell ( 1980[ ).) The dynamics of 11 are 
crucial for understanding the sensitivity and predictability of the underlying 
model. 

In his opus magnum, which was published in Ukraine in 1892, Lyapunov 
presented his theory of stability for finite-dimensional ordinary differential 
equations. This work includes his study of the non-autonomous linear prob- 
lem (10), see Lyapunov (1992) (yes, 100 years later.) One of Lyapunov's 



goals was to develop an analogue of the well-known eigenvalue-eigenvector 
theory, for the solutions of the autonomous problem, to the study of solutions 



of general non-autonomous equation (10). 



The approach developed by Lyapunov begins with the 4 Lyapunov Rela- 
tions of exponential growth: 



limsup - log(||f/(6',t)xo||) and 

t— ^iboo t 



liminf-log(||?7(^,t)xo| 

t— ^itcxD t 



where Xq ^ 0. Lyapunov was interested in, as are we, the case where these 
four limits are equal, and 



A(^,a;o) =^ lim ^ log(||f/(^^, t)xo||) 

t— s>— oo t 



lim^log(||f/(e,t)xo||), (14) 



where xq ^ 0. The linearity of U{9,t)xo implies that s\{9,xo) = \{9,sxo), 
for s 7^ 0. Hence one can assume, as we do, that xq is a unit vector, i.e., 
Ilxoll = 1. When (14) holds, then X{9,Xq) is a strong Lyapunov exponent, 
and the unit vector xq is an (associate) Lyapunov vector. The Lyapunov 
spectrum, LYE, is the collection of all such A(^, a;o), with G 21 and ||a;o|| = 1- 
For example, if A{9 ■ t) = Aq is an autonomous matrix, then the Lyapunov 
spectrum consists of all real numbers A that satisfy A = Re z/, where u is an 
eigenvalue of Aq. 
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A Lyapunov vector is not an isolated vector, rather it spawns a line of 
Lyapunov vectors (through the origin) in R-^. That is to say, a Lyapunov 
vector is a point in V^~^, the {K — l)-dimensional projective space. For a 
given vector v ^ in M^, we will use [v] G V^~^ to denote the unique line 
in that contains v. (Note that [v] = [—v].) Conversely, when one maps a 
line [u] in V^'^ to a vector u G M'^, we require that the pre-image u lie on 
the line [u] and that \\u\\ = 1. One should note that V^"^ is a metric space. 
The projective metric dpj-{ui,U2) is defined for nonzero vectors and u"^ in 
by 

dpr{u^,u^) = min llsin"*^ ± S2M^||, (15) 

Sl,S2 

where Si and S2 are real numbers that satisfy = ||s2M^|| = 1. This 

metric is used for measuring the distance between the lines [u^] and [u"^] in 

-pK-l 



Since the solution operator U{6,t) of the linear problem (10) maps lines 
in onto lines, one can use this operator to define a related projective flow 
Ti{6,t) on V^^^^ by means of the relation 

S(e,t) [u] = [U{e,t) u], for [u] G V^-\ 

Using this, one obtains an equivalent fiow [11] on x 21, where 

[U]{t){9,[xo]) = {e-t,J:{e,t)[xo]), hi {9,[xo]) e^xV""-', (16) 



compare with (12). One obtains additional information about the dynamics 
on the projective fiow [n(it)], by using the Lyapunov vectors, as is noted 
below. 

2.4 The Finite-Time Lyapunov Vectors 

Next we turn our attention to the question of finding good finite-time ap- 
proximations of these Lyapunov vectors. 

We begin by constructing a piecewise autonomous approximation of the 
Tangent Linear Equation for ([T]). To this end, we replace (10) by 



dtY{t) = An Y{t), for tn<t< t„+i, (17a) 

r(0) = Fo, (17b) 
tn = n6t, n = 0,l,2, ■■■ ,iV-l. (17c) 



where A„ = A{y{tn)) := j^\y={y{tn)) y{tn) is the value of an exact solution 
of ([T]) at t = The solution of this system at the grid points tn is, explicitly 
and recursively, given by 

= e-^*^" F,, forn = 0, 1, 2, ■ ■ ■ , iV - 1. (18) 
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Consequently, Y/v, the solution at time T = N 6t, is given by 



Y{T) = Yn = Z{T)Yo, (19) 

where 

We now define the approximation of the Lyapunov Vector associated with the 
largest Lyapunov Exponent Ai, at the time t = T - the finite-time Lyapunov 
vector, as the direction of steepest ascent for the matrix Z{T). For example, 
if one had used an explicit Euler scheme, then W{T) would be an Eulerian 
approximation of Z{T), where 

W{T) = {I + 6t An^i) ■■■{I + 6tAi) ■■■{! + 6tAo). 

The finite-time Lyapunov Vector is the singular vector corresponding to the 
largest singular value of Z(T). 



In practice (see e.g. Section 4.2 for the case of the CY92 model), the 
finite-time LV will be computed by directly solving a discrete approximation 
of the LTM and rescaling the output (the rescaling can be done at arbitrary 
intervals of time, since the problem is linear.) 



2.5 The Search for the Red Vector 

One of the main contributions found in the JPS87 manuscript is an indepth 
study of the interactions between two major theories of the longtime dy- 
namics of nonautonomous, linear differential systems. Dynamics of nonau- 
tonomous, linear differential equations: 

• Exponential Dichotomies (and Continuous Foliations) and 

• the MET (Multiplicative Ergodic Theorem) and Ergodic Measures. 

We view the JPS87 manuscript as a toolkit to be used in the analysis of 
the dynamics of related linear systems: It is this united theory, as we now 
show, that forms the mathematical foundations of the theory and applications 
of bred vectors. We begin with the first aspect: Exponential Dichotomies and 
Continuous Foliations. 

Consider the family of shifted semifiows 

Ux{e, t) ^= e"^* U{e, t), for A e M, 
and the associate skew-product fiows 

Ux{t){e,xo) =^ {e-t,Ux{e,t)xo). 



13 



Let Ai he a compact invariant set in 21. As noted in Sacker and Sell (1978), 
the skew-product flow 11^ is said to have an exponential dichotomy over Ji4, 
provided that there exist projectors Pa and Qx and constants Kq > 1 and 
a > 0, such that Px{0) + Q\{6) = I, and such that, for all 6' G and 
u G M^, one has 

\\Ux{0,t)Qx{e)u\\ < 7^0 e-"* II u II, t > 
\\Ux{0,t)Px{e)u\\ < Koc'^'Wul t<0. 



When there is an exponential dichotomy and (2.5) is valid, then the stable 
and unstable linear spaces, Sx and Ux - which are respectively the ranges 
of the projectors Qx and Pa - satisfy important dynamical properties that 
describe the exponential growth rate of selected solutions. For example, the 



system (2.5) is equivalent to: 

\\Ux{0,t)u\\ < Koe-'^'Wul for t > 0, n G n{Qx{e)), 
\\Uxi0,t)u\\ < Koc'^'Wul for t < 0, M G 7^(PA(^))• 

Let SSS denote the SS Spectrum (aka Sacker-Sell Spectrum) for 11, which 
is defined as the collection of all A G M such that IIa does not have an 
exponential dichotomy over Ai. 



The Spectral Theorem in Sacker and Sell (1978) describes the continuous 



foliation, and other properties of the flow 11 over Ai. More precisely, there 
is an integer i, where I < i < K, such that SSS is the union of £ closed, 
bounded intervals, that is. 



e-i 



SSE = y [ai_i,bi_i], with a£_i < bi^i < a^. 

1=0 

Also there is a continuous foliation 



■i-l- 



i=0 



where {Ve{9), ■ ■ ■ , Vi(6')} is a linearly independent, continuous family of sub- 
spaces of with Ei=o dim Ve-i{e) = K,hT 6 e M. As is shown below, 
the right-most interval [ai, 6i] plays a special role in the study of bred vectors. 

For 1 < k < i, we let [a^, bk] denote the /c*'^ spectral interval and let 
Vk{9) denote the corresponding subspace given by the continuous foliation. 
We next explore the important connections between the exponential growth 
rates of solutions with initial conditions in Vk{6) and the k^^ interval [ak, bk\- 
Among other things, we will encounter the Monotonicity Property and the 
Strictly Monotone Property. 
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Let V and A be real numbers that satisfy z/ < A, and both U,y{6,t) and 
Ux{6, t) have exponential dichotomies over Ai. Then the Monotonicity Prop- 
erty holds: 

C Sx, while Ux C U^. (21) 

The Strictly Monotone Property is a consequence of the observation that the 
following three statements are equivalent: 

• One has S^, ^ Sx- 

• One has Ux ^Uv- 

• There is a spectral interval [a^, bk] in the interval (z/, A). 

This brings us to a basic property. Let be+i and qq satisfy: bg+i < ae and 
61 < oq. Next we fix u and A so that 



bk+i < u < ak <bk < X < ttk-i. 

Then neither u nor A lie in the Sacker-Sell spectrum SSTj. Furthermore, the 
interval [u, A) contains the spectral interval [ak, bk]- By the Strictly Monotone 
Property, the space Sx is larger than Si,, while Ui, is larger than Ux- Moreover, 
as is shown in Sacker and Sell (1978), one has: 



u,{e) n Sx{e), 



for all eeM- 



(22) 



It is a consequence of the relations (|22j) and (20) that if the initial condi- 
tion {6,Xq) satisfies Xq G Vk{6), then the solution U{6,t)xQ satisfies 



\\U{d,t)x4 <Koe 
\\U{e,t)xo\\ <Koe 



iX-a)t\ 
{u-a)t\ 



Xo\ 



for t > 0, 
for t < 0. 



(Note that Kq depends on the choice of u and A.) 

It should be noted that all the terms used above, including Pa (6*) and 
Qx{0), vary continuously in 6- Furthermore, the exponential dichotomy is 
robust, in the sense that it varies continuously under small perturbations. 
Small changes in the model result in a related exponential dichotomy with 
small changes in P^, Qa, Kq, and a, see Pliss and Sell (1999). 

Moreover, it is shown in Sacker and Sell (1976) that if A e SSS, then 
there is a {0,xo) G x IR-^, such that ||a;o|| 7^ and the solution Ux{9,t)xo 
satisfies: 



sup \\Uxi 



, t)xo|| < 00. 



(23) 



What will become apparent shortly is that the red vector must be in the space 
Vi{9). Furthermore, since A = 61 is in SSS, the pair {9,Xo), that arises in 
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(23) for this choice of A, is a candidate for the "red vector" designation. A 



red vector must satisfy (23), but the converse need not be true. More on this 
later. 

As noted above, the second tool to be used in the theory of Lyapunov 
exponents/vectors is the Multiplicative Ergodic Theorem (MET) and the 
ergodic measures on Ai. One finds in JPS87 a study of the links between 
the LYE and the SSS. While the SSE leads to a continuous foliation, as 
noted above, the MET leads to a "measurable" refinement of this continuous 



splitting, as is noted in Remark 4.2.9 on pages 177-179 in Arnold (1998). 
The latter reference is noteworthy because it contains various extensions of 
the MET to problems not originally envisioned in the pioneering works of an 
earlier generation. 

The Multiplicative Ergodic Theorem: Let be a (non-empty) compact, 
invariant set on the attractor 21, and let n be an ergodic measure on Ai with 
/i(A^) = 1. Then there is an invariant set Ai^ in A^, with ii{M.^) = 1, and 
there is a k, with i < k < K, such that the following hold: 



1. There is a measurable foliation 

j=k^i 

= Wk-,{e), for 9 e M^, 

j=0 

where {Wk-j{0), ■ ■ ■ , Wi{9)} is a linearly independent, measurable fam- 
ily of subspaces of with dim Wj{6) = rrij > 1, for 1 < j < 
k, (mi + • • • + mfc) = K, and all 9 G A^^. 

2. There are real numbers A^, for 1 < j < k, that are the strong Lyapunov 
exponents A(^,Xo) = Xj, for all {9,xq) G Ai^ x Wj{9), with ||a;o|| = 1, 
and one has A^ < ■ ■ ■ < Ai. 

3. The Ergodic Spectrum LYE(/i), which depends on the ergodic measure 
11, is this collection {Ai, ■ ■ ■ A^}. The Lyapunov Spectrum LYE is the 
union 

LYE = U^LYE(^), 
over all ergodic measures /x, is used below. 

4. For each j, the measurable vector bundle 

M, X Wj{9) = {{9,xo) eM,x W,{9)} 

is an invariant set for the projective flow. Furthermore, because of 
the exponential separation between these vector bundles, the bundle 
Ai^ X Wj{9), with j = 1, is an attractor for the projective flow. 
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Connections between the SSS and LYE: As is shown in JPS87, the fol- 
lowing relations hold: 

• One has LYS C SSS. 

• For each i with 1 < i < i, there is an ergodic measure fi with support 
in the spectral interval [ai,bi]. 

• Assume that the invariant set Ai is "dynamically connected", that is, 
Ai cannot be written as the union of two disjoint, nonempty, closed, 
invariant sets. Then the following holds: For bi, the largest value in 
SSE, there is an ergodic measure /i with the property that A = 61 is 
a strong Lyapunov exponent. It follows that Ai = 61 is a Lyapunov 
exponent for some ergodic measure /i on A^, and consequently there is 
a unit vector xi in Wi{6) with the property that xi is a red vector. 

• The previous item is valid for any of the endpoints {oi, ■ ■ ■ , a^; 61, ■ ■ ■ bk}, 
but the related ergodic measures may differ. 

It can happen, as in the case with the CY92 model, that there is a unique 
red vector. Furthermore, in this CcLSG, cLS IS noted in Section Hi one has 
dim Vi{9) = 1. Hence one has Wi{9) = Vi(^), and there is a unique red 
vector in the projective flow. Moreover, due to the exponential dichotomies 
occurring in the CY92 model, the red vector is robust, and it varies contin- 
uously with small changes in the model. 

On the other hand, when dim Wi{6) = 1, there is always a unique red 
vector. If in addition, one has dim Vi(^) > 2, then the red vector may be 
only measurable and not continuous. In short, the red vector need not be 
robust. 

3 Continuum Limits of the BV and EBV Al- 
gorithms 

We now elucidate further the relationships between the BV and EBV algo- 
rithms and the dynamics of the underlying system. In particular, we now 
formalize the connections between these algorithms and the solutions of the 
linear tangent equation, i.e., the finite-time Lyapunov vectors. We are in- 
terested in the behavior as the step size 6t goes to 0. To accomplish these 
aims we revert to a continuum formulation and for simplicity, assume that 
the vector y{t) in ([T]) is a member of the Euclidean space M"^. Since energy 
norms are used in many geophysical fluid mechanics problems, we will take 
the usual /^-norm. We denote the norm and inner product by || ■ || and (■,■), 
respectively, with || ■ |p = (-, ■) 
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The notation in this chapter differs from what has been set in the rest of 
the paper in minor ways, hke the letters representing the functions. This is so 
to emphasize that unhke what one computes in practice, the dynamical sys- 
tems here are continuous in time. However, everything is clearly explained to 
avoid ambiguities without cluttering the presentations with technical details. 

We stress that the term continuum limit refers to the rescaling time 
(sometimes called a cycle), not to the numerical integration time step. In 
a numerical context, this corresponds to a strategy where both integration 
time steps and the rescaling times are small. While this has no importance 
for linear systems, it leads to different outcomes when applied to nonlinear 
systems even for quite small perturbations amplitudes. 

We first recall the system Q: 

dy 

-dt=^^y^^ (24) 

1/(0) = yo. 

Our first goal below is to obtain the formula 

W{T- e) = W{t^- e) + f [G^ {W{t- e),t) - (G^ {W{t- e),t), W{t; e)) W{t; e)] dt, 

Jto 

(25) 

where for all t >to, 

W{to; e) = Wo; \\Wo\\ = 1; {W{t; e),t) = ^Gy {eW{t; e),t), 

and Gy{6y,t) = G{y + 6y) — G{y) for the limiting case 6t ^ 0, to < t < T. 
The vector W(t;e) corresponds to a continuously rescaled bred vector V{t) 
at amplitude e and with initial perturbation V{to) = eWo- {W has the same 
direction as V, but it has amplitude 1). We will streamline the presentation 
by skipping some of the techical details in derivation, and we will assume 
at the outset that G in (24) has the necessary differentiability properties to 
make all the mathematical steps rigorous. 



Let y and y + 5yhe two solutions of (24) that satisfy the initial conditions 
yito) = yo and {y + Sy){to) = yo + ^yo- Then 6y(t) is a solution of 

j^6y = [G{y + 6y) - G{y)] , 6y{to) = 6yo. (26) 



Even if G is autonomous, the resulting equation ( 26 ) for 6y is non-autonomous. 
Integrating (26) for t > to, we obtain 

r-to+St 



/ rio+dt \ 

Sy{to + 5t) = (^[y{to) + 5y{to)] + G{y + 5y) dtj 



to+St 

y{to) + / G{y) dt 

to 



(27) 
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We next define e = \\Sy{to)\\, and a family of vectors v{t) on a sphere centered 
at zero with radius e in M.^ by 



Assuming that 6y is bounded away from 0, v{t) obeys the evolution equa- 
tion: 



dv 



\\6ym dt \\6y{tW 
1 



GyiSy,t)\\Syit)\\ 



1 6y{t) d 



(2J 



2\\6y{t)\\ dt' 

which makes explicit the fundamental dependence of the time evolution of v 



on the norm. Taking the scalar product of (26) with 6y{t), we also have 

1 d 



2dt 



\\6yr = {Gyi6y,t),5y{t)). 



Substituting this relation in (28) gives after some simplifications, 

1 dv ^ Gyi6y,t) _ ( Gy{5y,t) v{t) \ v{t) 

mmdt \\6ym \ \\6ym 'mm) ¥ym 



(29) 



We integrate ([29j) independently on successive intervals 

[to, to + 5t), [to + 5t, to + 25t), . . . , [to + (AT - l)5t, T) 



where T = to + N5t and denote the solution of (29) on the interval Ik 



[to + (/c — 1) (5t,to + kdt) by ffc(t). This is a sytem of N integral equations 
with free parameters (the integration constants) f^, /c = 1, . . . , A^: 



{k-i)5t \ \\'jyk{T')\\ 
to + (A; - l)6t < t < to 



Gy{6yk,T)) f Gy{ 6yk,T) Vk{r) \ Vkir) 
kiT)\\ ' H\\ 
kSt, k = l A^, 



c/r, 



(30) 



where Syk satisfies the (27) with the initial condition at t = to + {k — l)6t, 
i.e., the left-hand boundary of the interval 7^. We observe that, at t = 
to + {k — l)6t, Vkit) is an approximation of the discrete bred vector with 
short rescaling time step 5t. Hence, we take it as the basic approximation 
for the continuously rescaled bred vector we seek. Note that the integrand 

Gy{&y,t) 



in (|30|) is simply the component of the vector jf/^i'^m*! perpendicular to ffc(t) 



m 



at time t. 
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We take 



i;° = 5y(to), vl= lim i;fc„i(t), l<k<N, (31) 

assuming the limits exists, and extend the domain of definition of Vk by 
setting Vk{t) = if t ^ 7^. Similarly, we extend the domain of definition of 
5yk. We then let Vn = vi + ■ ■ ■ + vn, and note that Vn is continuous in 
time on [to, to + N5t\ (extended to t = to + N5t by continuity) and piecewise 
differentiable. We define 5Y^ analogously as the sum of extended 5?/fc's. Then 

V = lim Vn = lim 

N^oo (5t-s>0 

exists and is continuous in time, for instance, when G has the necessary 
differentiability properties so that the sequence Vjy, possibly after passing to 
a subsequence, converges uniformly on compact time intervals, as iV — )■ oo 
(or, equivalently, as 6t — )■ 0). Since V^v coincides with the right-hand side 
limit of 6Yn at every grid point, the difference between Vn and 6Yn goes to 
zero as iV — )■ cxD in an appropriate norm that is allowed by how smooth G is. 
When G is continuously differentiable, for example, this convergence would 
be uniform on compact time intervals. Thus, by passing to the limit N ^ oo 



in (30), we obtain the following representation formula for V: 



V ^ VM + £ [oAVm) - {G,(V(t),t), ^) ^) M. (32) 

Under the regularity assumptions above on G, the integrand is continuous 
and hence is a mild solution of an associated differential equations. There- 
fore, we can bootstrap and prove the further regularity of V. Note that, by 
construction, 

e=\\Vito)\\ = \\V{t)l 



Therefore, Equation (32) is equivalent to (25), if we take 6yo = eWq. 



An immediate consequence of (25) is the following. When G is sufficiently 



regular, W{t;e) converges uniformly on [to,T] as e — )■ 0. Let us denote the 



pointwise limit by W{t). Then, by (25), and the observation 



1„„ Gl (W(t), t) = Mm Gfa + ^H--W)-Gfa) = z)„,„G(,). 

where the right-hand-side is the directional derivative of G in the W{t)- 
direction, it follows that W{t) satisfies the differential equation 

^ = A{t)W - {A{t)W, W)W, (33a) 
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where 
and 



A{t) = DyG\ 



y^\y=y{t)} 



We recall the Linear Tangent Equation for the problem (24): 



~dt 



A{t)X, A{t) = DyG\ 



y^\y=y{t), 



(33b) 



(34a) 



X{to) = Wo, X eW\ (34b) 

where DyG is the Jacobian matrix of the map G{y), and y{t) is the solution 
of (24). As in Section |2| U{t, s) denotes the solution operator of (34), taking 
the solution at time s to the solution and time t. With a slight abuse of 
notation, we write U{t) := U{t,tQ), so that in particular the solution X{t) of 



(34) at time t is simply U{t)Wo. We recall also that U is a semiflow, i.e., 



U{t + s) = U{t + s,t)U{t), Vt, sgM. 



(35) 



What is the relation between (33) and (34)? It can be easily seen that 



rescaling X{t) to unit length and restarting the integration of (34) with 
X{t)/\\X{t)\\ at any t„ € [to,T] does not alter value oi X{t)/\\X{t)\\ for 
t > tfn- By induction, this follows for any finite number of rescaling-restarting 



cycles. This is because of the linearity of (34). We leave the details to the 



reader. By a direct approximation argument, or by a similar argument that 



led to ( 33 ) , we obtain 



\m\\ 



Wit), te[t,,T]. 



Thus, we have just established that, the continuously rescaled bred vector 
W(t;e), i.e., the solution of (25), converges to the corresponding solution 



(rescaled to size 1) of the tangent linear equation (34) with the same initial 
data, i.e., to the solution of (33), as e — )■ uniformly on the compact set 



[to, 7"], assuming that G is smooth enough. Since the solutions of (33) have 



constant magnitude 1, it is the linear analog of (25) in terms of continuous 
rescaling. 

While continuous rescaling is inconsequential for linear equations, this is 
not the case for the nonlinear ones. However, when the initial amplitude is 
very small, the results above indicate that one might be able to rescale less 
often and still get comparable results. This is simply because it takes longer 
for nonlinear effects to start to dominate the picture. 
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3.1 Continuum Limit for EBV 



We discuss the continuum limit for the EBV algorithm briefly, by confining 
ourselves to a short, heuristic description. Writing a system of differential 
equation for the EBV algorithm is not as straightforward as for the BV 
algorithm. The perturbations Sy{L) in the ensemble that grow fastest at the 
largest amplitude (before rescaling) at time t will obey the equation (25) as 
long as they are the top contenders. All the other perturbations will satisfy 
the following differential equation: 

^^=G,(MO,t)-— (36) 

where l El, and dh{t)/dt is the maximum growth rate of the norm among 
all ensemble members with the maximum amplitude e at time t. Assuming 
enough smoothness on h{t), (36) can be heuristically obtained by taking 

the right time derivative of the vector 6y{L){t)j-^, where h{T),T > t is 

the maximum norm of the ensemble members at time r in absence of any 
rescaling after time t (If Equation 26 were in charge). As e approaches zero 
for time t fixed, all ensemble members point in directions along which (36) 
approaches the linear tangent equation, since dh{t)/dt — )■ by continuity 
of Gy at zero with 6*^(0, t) = 0. A detailed treatment complete with the 
technical aspects will be presented elsewhere. 



3.2 Tuning the Maximum Amplitude 

We note that the convergence to the linear tangent equation is expected 
to be quicker for the EBV compared to the BV for most of the ensemble 
members. This is manifested in the equations, and later verified numerically 
in the next chapter. The acceleration is due to the inherent size ordering in 
the EBV algorithm: information on the scales where linear tangent equation 
dominates tends to be preserved in a robust manner against the changes in 
the parameter e. On the other hand, at finite amplitude away from 0, the 
rescaling rule in the BV algorithm might create or prolong the life of certain 
instabilities due to the nonlinear effects. The algorithm output needs to be 
examined independently to check whether these vectors are in fact relevant to 
the dynamics or just artefacts of the rescaling strategy. We believe that the 
EBV algorithm is more resistant and robust in this regard. The parameter 
e is tuned under different considerations in the EBV and BV algorithms. 



3.3 Separation of Scales 

In nonlinear systems, it is possible that a perturbation will not grow very 
rapidly in the zone of perturbations with size near e, but instead, a multiple 
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of this perturbation can be dominant among the smaller perturbations and 
might grow quite fast there. It is also possible that the dominant vectors 
in different amplitude zones will be close in structure, but very different in 
growth characteristics. These theoretical considerations are realized in the 
Lorenz63 system quite strikingly. See Figure |4] and its discussion, especially 
the recurring patterns of "spears" therein. This can be seen very easily from 
the following relation: 

\\Vitn)\\ ^ ll^(^n)|| 

where R™° is defined similarly to ([9]), V{t) denotes a EBV(t) member, and 
V{t) denotes the same vector right before the rescaling. In a nonlinear sys- 
tem, this ratio can be larger than 1, for some range of perturbation sizes 
smaller than e. This is the main mechanism that creates the vector zones 
of different magnitude with zonal growth characteristics, or a separation of 
scales. 

Another case in point will be presented in Section |4] for the CY92 model. 
Even at late times and small amplitudes, we can see perturbations surviving 
the rescaling strategy and they resemble what one would get from the usual 
BV algorithm at those sizes (and the finite-time Lyapunov vectors) very 
closely, whereas the dominating perturbation of size e resembles a particular 
BV of size e. This example actually illustrates more. There are members of 
the ensemble in small magnitudes that are slow in aligning with the dominant 
directions. 



4 Applications of the Bred Vector Algorithms 

We will be comparing the BV and the EBV on two problems: The Lorenz 



equations, or Lorenz63 (see Lorenz (1963)), and a dissipative and forced non- 



linear partial differential equation that arises in modeling the thermohaline 
circulation (see Cessi and Young (1992)). The latter will be denoted as the 



CY92. It is a Cahn-Hilliard equation and it will shown to have an inertial 
manifold. We will also have occasion to compare the BV and EBV results 
to the finite-time Lyapunov vector outcomes. 

Throughout we will use an explicit fourth order Runge-Kutta time march- 
ing scheme, for the calculation of the base solutions as well as for the calcu- 
lations of the BV, EBV, and finite-time Lyapunov vectors. 



4.1 The Lorenz63 Model 

The finite-dimensional, nonlinear Lorenz63 model has often been used as 
benchmark for testing sensitivity and, in particular, as a test problem for 
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BV (see for example Evans et al. (2004)) 



Let X = (Xi, X2, X3) G and let 

A=l r -10 , N{X) = -X1X3 , DNiX) = ^iV(X). 
\ —b J \ X1X2 J 

(37) 

DN{X) is a 3 X 3 matrix- valued function. The Lorenz model is described by 
the solutions of the nonlinear system: 

dtX = AX + N{X). (38) 

The associated Tangent Linear Model is the skew product system 

d,X = AX + N{X) 

dtU = iA + DN{X)) U, ^ ' 



where U = (?7i, f/2, f/3) G The [/-equation in (39), which is a linear 
equation, is of special interest to us. 

We will use the notation for the Tangent Linear Model introduced in Sec- 
tion |2^ With the initial condition 9 = Xq e ]R^ we let 61 - t = S{9, t) = X{t) 
denote the solution of the nonlinear equations (38) that satisfies S{6, 0) = 6. 
In this study we set r = 28, 6 = 8/3, a = 10. It is well-known that for these 
values of the parameters the Lorenz63 model has a chaotic global attractor. 
For the non-autonomous [/-equation; dtU = {A + DN{6 ■ t)) U, we let 

U(t) = U{9,t)Uo, denote the solution operator, (40) 

where f/o G M^ and f/(0) = Uq. 

We will use the term "attractor" to refer to a compact, invariant set 21 
that attracts a neighborhood of itself. As a result, 21 is Lyapunov stable, as 
a set. Consequently, for < e < 1, there is a family of £- neighbor hoods, 
iVe, of 21, where each neighborhood is positively invariant and 21 = He^oNs. 
When we write that 6 is near 21, we mean that 9 & N^, for some small e > 0. 



See Chapter 2 in Sell and You ( 2002 ) , for a history of this concept and more 
information. 

By using a somewhat different -but equivalent formulation- Toth & 
Kalnay have suggested that the time evolution BV{t) is a good approxi- 
mation of the tangent linear solution f/(t), over bounded time-intervals. As 
a consequence of the continuum limit theory in Section [3} we see that this 
perceptive observation has a solid mathematical basis, provided that the 
time-step 5t is small and the perturbation amplitude is small, as well. 

In order to find a numerical validation of the continuum limit theory 
described in Section [3| we will calculate the distance 

D{t,5t) := dpr([W(t)], [f/(t)]), (41) 
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where dp^ is the projective metric on 21 x V^'^ (see (15)), and BV{t) is a BV 
at time t. For this calculation, we assume that both BV{t) and U{t) satisfy 
the same initial condition (yo, 53^o) at t = 0, and that uq = 9 is near the 
attractor 21. We also use identical ensembles of perturbation vectors for both 
the BV and the EBV algorithms. We then fix T > 0, and we examine the 
distances D(T,6t) (in the projective space), for different choices of 6t. Our 
goal is to show that D{T,6t) becomes smaller, as 6t gets smaller. The max 
and min values, for both the BV and the EBV algorithms, are reported in 
Table [1} For the calculations used to generate the data in Table [T| we had set 
the initial conditions for both BV{t) and U{t), at t = 0, so that X = 0.5688, 
Y = 0.4694, Z = 0.0119, and 6yo = [1, 1,0]. For the ensemble, we took the 
3-fold Cartesian product of the set 

{-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1} 

in M}, and projected this product onto the sphere of radius 0.1, centered 
at the origin, in M^. By eliminating the repetitions introduced with this 
projection, one obtains an ensemble of 584 distinct points in M^. (Note that 
the initial conditions are "near" the Lorenz attractor.) The trajectory for the 
nonlinear problem was computed with a step size dt = 0.0001, so that any 
observed variations are only due to the differences in the algorithms and the 
relative effects on the perturbations. In Table [TJ we fixed T = 2 and made 
two choices for 5t, namely 6t = 0.004 and 0.001. The results compare the 
outcomes of the BV and the EBV algorithms. As noted in Section 3, both 
algorithms EBV{t) and BV{t) approximate solutions of the tangent linear 
equation, as 6t goes to 0. By using the projective metric one can compare 
the rates of convergence in terms of this metric. 

The maximal values of D{T,6t), for both the BV and EBV algorithms, 
are essentially the same for both choices of St, and the minimal values are 
essentially the same for 6t = 0.004. When one moves from 6t = 0.004 to 
6t = 0.001, both minimal values decrease, as expected. However the drop 
in the minimal value for the EBV algorithm is substantially larger than the 
drop for the BV algorithm. By using the perturbation corresponding to the 
minimal drop in D{T,6t), for each of the algorithms, one arrives at the best 
approximation given by for the given algorithm. Clearly the min BV and 
min EBV columns in Table [T] shows that the EBV algorithm yields a better 
approximation than the BV algorithm. 

One can, of course, use different initial conditions for the two solutions 
BV{t) and U{t) and/or larger perturbations. However, one cannot expect to 
replicate the results seen in Table [T] in that case. First, there is a transient 
phase, which ends when the two solutions are "near" the attractor. Even if 
this transient phase is short, one still has a problem. While the attractor 
is stable, as a set, one still has a problem because the flow on the attractor 
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Table 1: Maximum and minimum values of D{T, 6t) with T = 2, for two 6t. 
See gip. 



6t 


max BV 


mill BV 


max EBV 


mill EBV 


0.004 
0.001 


1.40 X 10-1 
7.77 X 10-2 


5.26 X 10-^ 
1.03 X 10-^ 


1.01 X 10-1 
3.48 X 10-2 


1.10 X 10-3 
4.16 X 10-5 



is generally not Lyapunov stable. Essentially all pairs of nearby orbits may 
diverge over long time intervals. 

An interesting outcome of the use of the EBV algorithm on Lorenz63 is 
that it exposes the fractal behavior of the Lorenz attractor. We computed 
a 584-member EBV, using the same parameter values, initial conditions and 
base trajectory as was used in generating Table [T| for the 6t = 0.001 case. 
We constructed plots by merging all the vectors between time 24 and 30, at 
intervals of dt = 0.1. In that time, at any given instant, one only observes 
a couple of members reaching highest amplitude. By rotating the same plot 
about the = Z axis, one obtains three views of the EBV's shown in 
Figure [1} Zooming into Figure [1^ by a factor of 8, 32, and 60, respectively, 
we obtain Figure [2] The fact that one observes similar patterns at several 
levels of magnification meets the requirement of 'fractal behavior' that was 
introduced by [Mandelbrot (1977). 

Figure [3] shows the time evolution of a sample BV and the corresponding 
finite-time Lyapunov vector with the same initial perturbation. The results 
were obtained with a time step of 6t = 0.005. The computations of the 
BVs and finite time Lyapunov vectors were performed with a time interval 
equal to the dynamics time step. At t = 0, X = 0.1493, Y = 6.2575, 
Z = 1.8407. The initial perturbation vector 53^o was [1,1,1]. As expected, 
small to moderate-sized perturbations produces similar results in the finite- 
time Lyapunov calculation and the Bred Vector calculation. 

Figure |4] is devoted to the results of a calculation of EBV{t) alone, for 
the same case. The figure depicts a rescaled time evolution of the norm 
of 98 distinct ensemble members in the EBV calculation. The rescaling is 
done with respect to the usual Euclidean t^-noim. The initial perturbations 
are made to sample a perturbation sphere of amplitude 1 about the initial 
conditions. The figure highlights the rapid decay of many of the vectors 
and the eventual size-ordering that is inherent in the EBV algorithm. For 
< t < 13, the results in Figure |4] describe the transient behavior and 
are rather chaotic. However, for t > 13 a very interesting pattern evolves: 
we see that the largest member of the ensemble takes on the value 1, for 
all t > 0, which is expected. (This corresponds to the Lyapunov vector 
with the largest exponent.) From the totality of the information on the 
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Figure 1: Viewed from three different angles, (a)- (a). EBVs, for times 24 
through 30, taken at 0.1 time intervals, for Lorenz63. Parameters and con- 
ditions are those used to generate Tahle\^ 
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Figure 2: Close-up views of Figure^^: (a) zoomed in 8 times; (b) zoomed in 
32 times; (c) zoomed in 60 times. Note the preservation of structure even as 
we zoom in. 
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Figure 3: The three components of BV and of the finite-time Lyapunov vec- 
tors, as a function of time, for the Lorenz63 system. (For parameters and 
initial conditions, see text). The vectors are nearly coincidental over the 
whole time span. The perturbation size is 1, in all of the components. For 
our choice of initial conditions the vectors exhibit more temporal regularity 
at earlier times. This transient behavior disappears att = 13, approximately. 

ensemble, it is possible to extract the next two Lyapunov vectors, where 
the exponents satisfy A3 < A2 < Ai. However, this would require a deeper 
analysis, and it is deferred to a future work. In fact, viewed as a whole, 
the graphs corresponding to those smaller than 1 have very useful structural 
information: Notice the recurrent 'spear-like' pattern, which occurs after 
t = 17. What is happening is that a group of 'small' vectors in the attractor 
grow rapidly, as they go around the horn in the Lorenz attractor. We propose 
that the rationale for this behavior is that the nonlinear equations of motion 
temporarily overwhelm the uniform rescaling rule for these small vectors. 
For t > 17, we are beyond the transient zone, and we would not expect such 
behavior to occur if the linear term AX strongly dominates the nonlinear 
term N{X) in the vicinity of the attractor. This is an excellent illustration 
of the fact that EBV algorithm preserves the role of the nonlinear terms in 
the equations of motion. 

This feature of the EBV can not be replicated by BV, even a BV with 
an ensemble of perturbations. The reason for this is that the rescaling rule 
for the BV algorithm forces all the perturbations to have the same norm, for 
each t > 0. Thus the Figure |4] for the EBV would be replaced by a figure for 
the BV, where all of the perturbations are plotted on the top line only. 
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Figure 4: Plot depicting the evolution of the 2-norm of an ensemble of 98 
EBV{t). Same parameters and initial conditions, as in Figure^ Initially, 
in the interval [0, 13] we note a very fast decay of some of the ensemble 
members, leading to a sorting in size, beyond that time. As described in the 
text, this outcome is one of the most distinguishing features of the EBVs, 
when compared to an ensemble of individual BV outcomes. 
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4.2 A Cahn-Hilliard Equation 



In their work on the thermohahne dynamics Cessi and Young (1992) pro- 



posed a coupled model for the circulation, salinity, temperature, and density 
of the oceans, with atmospheric forcing. The crux of the model is the par- 
tial differential equation for salinity: The slow-time dynamics of the ocean 
salinity S{x,t), zonally-averaged, and as a function of latitude x G [— 7r,7r] 
and time t > 0, is described by 



dS 

~dt 
S{x,0) 



a — [f{x) + fiS{S - sin(x))^ + S - 
So{x). 



t > 0, 



(42) 



The equation is subject to zero-flux and zero-stress boundary conditions at 
the poles, however, we will be considering periodic boundary conditions (for 
steady and periodic forcing /(x) as well as equilibrium solution sin(a;) the 
conclusions that follow apply to the zero-fiux case). The positive parameter a 
affects the strength of the linear stability of the model. We fix a = 3.5 x 10~^, 
7 = 0.001, and /i = a/10. The forcing f{x) is a prescribed function that 
reflects balances of evaporation and precipitation of freshwater; it can be 
symmetric, about the Equator [x = 0), but is more typically, asymmetric 
Eyink (2005)). The second derivative of the forcing function, will be 



see 



chosen to be 



dxxf{x) 



^ (36 — 39/i^ cos^ x — 81 cos^ x + 8yU^ + 25/i^ cos"^ x) sin x, 
for X G [— TT, 0], 



(3/i^ cos^ X — 3 — fi"^ 
for X G (0, vr]. 



(43) 



smx 



It is shown in Figure [5^. The initial condition chosen for this computation 
is 5*0 (x) = cos(x). The base solution obtained numerically is displayed in 
Figure |5|d. 

A great deal is known and can be said about the mathematical structure 
of CY92 and its solutions. By letting S = dxU and / = dxF the CY92 can 
be related to the Cahn-Hilliard Equation (CHE) with forcing: The general 
CHE, for u above, is 



dtu + v/y 



u 



A{g{u)) + F for X G and t > 0, 



(44) 



where u = u{t,x) is a scalar field, u > is a constant, A is the Laplacian 
operator, is the bi-harmonic operator, and g is a. polynomial of degree 3: 
di"^) — X]j=i with > 0. The term F = F{x) is the forcing function. 

In general, the domain VL may be an open bounded domain in the Euclidean 
space M"*, with 1 < m < 3. However, we restrict our attention, to the case 
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(a) 



(b) 



Figure 5: (a) Forcing function fxx, used in the CY92 model; the nu- 

merical solution to CY92 is shown in (h). See Figure^for the BV algorithm 
and finite-time Lyapunov vector algorithm outcomes. 

where m = 1 and fl is the interval Q = (— tt, it), with boundary F = {±7r}. 
As we will see, the multiplicity of the largest Lyapunov exponent for this 
problem is 1. 

For the analysis of the solutions of the CHE, one will use the standard 
Sobolev spaces = H^{Q), where j > is an integer and = L'^{Q). As 
usual, the inner product and norm on is denoted by (-, ■) = (-, ^ and 
II ■ II = II ■ II 0- For j > one uses 



\u\ 

'J-r^' II iij 

\a\=j+l ' 

where 

Qa 

and 

The finite-time Lyapunov vector algorithm for the CY92 model is derived, 
by first rewriting the equation as 

dtS = ad,, [/x^S (5 - r^f - rf{x) + S- ^^d,,S] =: F{S). 
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We opt here to first linearize and then discretize. To obtain the tangent 
hnear equation, we need to find the hnear map L such that 

F{S + bS) - F{S) = C{SS) + o{6S) as \\6S\\ 0, 

where || ■ || is a norm. Typically, this will be the norm in the space where 
solutions live, such as the Sobolev space for S, or dictated by physical 
considerations. In this infinite-dimensional model, different norms are not 
necessarily equivalent. 
First, we note that 

{S + SS} {{S + 6S} -^f-S{S- r^f = {S - r^) (35 -r)) + o{\\5S\\). 

We let L be a finite-difference approximation to C. The discrete tangent 
linear equation becomes the product of two matrices AB, where B is the 
matrix whose entries are found by discretizing the operator 

IJ?{S-'q) (3S-rj) + Id-^'d,,, 

where Id is the identity matrix, and A is the matrix corresponding to the 
discretization of ad^x- The finite-time Lyapunov algorithm is obtained by 
solving 

dtSe ^Lse^ABse^ A{B{5e)), 

for t e [0, T], subject to some initial vector perturbation 6Q{0) — s, where SO 
is the discretization of the perturbation SS. The perturbation is normalized 
to the norm of the initial perturbation at each 6t, the time step of the explicit 
Runge-Kutta 4 time integration scheme employed here. T = 70. In the 
simulations, 6t — 0.01. We applied second-order centered finite differences 
in space, and used 121 grid points, — tt = xo,xi, ...,xi2o — tt. 

By construction the finite-time Lyapunov vectors arc not amplitude sen- 
sitive, however the BVs are, thus different amplitude perturbations will yield 
different BVs, in the case of a general nonlinear problem. This outcome has 
important practical implications, if one would like to use BV to either infer 
the structure in the field, or the degree of sensitivity of the outcomes to per- 
turbations in initial conditions. A challenging problem could thus arise in 
the context of large-scale simulations: what is considered a large structural 
change or a highly sensitive outcome is physics-dependent, perhaps even dif- 
ficult to surmise quantitatively; the physics in question may not be fully 
understood and thus a reasonable perturbation amplitude is simply guessed. 

We ran the finite-time Lyapunov vector and the BV algorithms, using 
the same initial condition, forcing and perturbation. First, we examine the 
effect of the size of the amplitude of the perturbation on the outcomes: we do 
so by keeping the shape of the perturbation fixed, changing only the overall 
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amplitude. Figure [6^ and b are plots of the final BVs and finite-time Lya- 
punov vectors, corresponding to 2-norm 0.25 and 0.025 sized perturbations, 
respectively. When the perturbations are small the BVs and the finite-time 
Lyapunov vectors are qualitatively consistent and, nearly so, quantitatively. 
The outcomes shown here are typical of the general case, that is, the qualita- 
tive and quantitative disagreement grows with an increase in the amplitude 
of perturbations. 




Figure 6: Comparison of the finite time Lyapunov vectors and BVs. Effect 
of amplitude of perturbation on the CY92 Model, with initial conditions Yq = 
cos(a;). The initial perturbation was esin(x). Comparison of the finite-time 
Lyapunov vector and the BV outcomes at t = 70. (a) Corresponds to e = 
0.25; (b) to e = 0.025. 

The shape or spectral content of the perturbation mattered as well. The 
spectrum of the perturbation is clearly important when projecting onto spec- 
tral bases for a reduced representation. To illustrate this, we use the CY92 
model, with the same forcing as before and same initial condition. We exam- 
ine monochromatic sine wave perturbations with wavenumber j = 1,2, ...,6 
of the form 

5Yo = ej sm[jx + 3 (j - 1) exp(l)], (45) 

where all of the perturbation amplitudes remain the same, ej = 0.25. The 
choice of the phase is inconsequential: to show this we chose non-commensurate 



34 



phases among the sine wave components. In Figure [7^ we show the BVs as- 
sociated with each of the perturbations (j = 1,2..., 6), at t = 70. Figure [7|3 
shows the EBVs at t = 70. Even at t = 5, shown in Figure we aheady 
see an amphtude-ordered structure in the EBV. 

The space-time plot of all EBVs is plotted in Figure [7|1; we note the 
very short transient phase, lasting till about t = 4, followed by a structure 
which is clearly dominated by the largest member of the EBV. The EBV 
calculation was run with identical parameters to those used in Figure [7^, 



with the ensemble consisting of the same initial perturbations in (45). In the 
BV case we are getting outcomes that do not have the reductive appearance of 
the EBV. On the other hand, there is no ambiguity to the prevailing ensemble 
member in the EBV case. This ensemble member has a clear correspondence 
in structure to the path and its perturbation field. The BVs should eventually 
agree with the EBVs, nevertheless. It is safe to assume that this will happen 
only after a very long time, longer than the time interval that might be 
suggested by the base solution. 



We also compared BV, for each of the perturbations in (45), to the out- 
comes of the finite-time Lyapunov vector calculation. Figure [8] compares the 
BVs and finite-time Lyapunov vectors for each j wavenumber perturbation. 
The amplitudes are set to ej = 0.25, for j = 1,2, ...,6. The structure of 
the finite-time Lyapunov vectors, to within a sign, is qualitatively the same, 
regardless of the wavenumber of the perturbation. The BVs are not. If the 
perturbation was made considerably smaller the differences between the BVs 
and the finite-time Lyapunov vectors would become small, as expected: As 
shown in Section |3} the BVs and the finite-time Lyapunov vectors must be 
similar to each other, provided the perturbations are small enough. For larger 
amplitudes, the BVs might show more structure since nonlinear effects could 
play a role in the structure of the BVs. Apparently, e = 0.25 is already in 
the range of large perturbations of the CY92 about the solution chosen. It 
was not clear whether a chosen perturbation is large or small, based solely 
on the CY92 model itself: it was only clear to us after a comparison of the 
outcomes of the finite-time Lyapunov case and the BV case. 

4.3 Cahn-Hilliard Dynamics 

The Bi-harmonic Operator Bu = A'^u on Q and the linear Bi-harmonic 
equation: 

dtu + z/A^ u = (46) 

play a basic role in the study of the CHE. We assume for now that the 
Operator B satisfies either the non-flux boundary conditions: 

du , d(Au) , , 

— = and = 0, on F, 47 

an on 
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Figure 7: (a) At t = 70, the superposition of 6 BV simulations, the j*^ one 



corresponding to a sine wave perturbation, as per (45). (h) Cross section of 
the EBVs, at t = 70, run with an ensemble of perturbations (see 45)- We 
note that, unlike the BV case in (a), the EBV has settled into a structure 
that clearly reflects the dominant EBV ensemble member, (c) Cross section 
of the same EBVs, at an earlier time: at t = 5. Even for short times, the 
EBV already shows an amplitude- ordering of its vectors, (d) Superposition 
of all EBVs, as a function of time and space. 
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Figure 8: Wavenumber dependence of BVs and finite-time Lyapunov vectors. 
CY92 Model, cosine initial condition and non- symmetric forcing (see Figure 
Comparison of the finite-time Lyapunov vectors and the BVs for sine 
wave perturbations of size tj = 0.25. The perturbations are the individual 
sine waves in (45), (a) j = 1, (b) j = 2, (f) j = 6. 
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or the periodic boundary conditions, see Sell and You (2002). In the sequel, 
we will assume that the forcing function F in ( 44 ) satisfies F E H^, for some 
j > 0. That is to say, 

Fdx = 0. 



Thus, by integrating ( |44| , one observes that any solution u = u{t, x) satisfies 

u{t,x) = u(t), 



uo{x) dx 



1 



for t > 0. 



(4J 



One seeks solutions of the CHE (44) in the Sobolev space V^. A mild 
solution u = u{t) of the initial value problem is given by 



u{t) 



-vB,{t-s) [a((7(m(s))) + F] ds. 



(49) 



We will denote the maximally defined solutions of (49) by u{t) = S{t)uo 



With uq G Vq, this mild solution is uniquely determined, with S{t)uQ G Vq, 
for < t < T(uo), where < T{uq) < oo. 

Because of ( [48| , we see that, for j > 2, the spaces 

H^^ ■= {(f) e ■.^ = 0} (50) 



are positively invariant spaces for the solutions of (44). Thus we will focus 



only on solutions S{t) uq that satisfy uq = 0. We also denote the collection 
of stationary solutions of the CHE by 

Q ■= {uq e : S{t)uo = uo for all t G M}, 

and Qo = QnH^ = QnVo\ 

4.3.1 The CHE With F = 0: 

The basic problem with forcing F = is of special interest. The equation 



(44) becomes: 



dtu + z/A^ u = A{g{u)) foixefl and t > 0. 



(51) 



To study the solutions of (|51|), one uses the Landau-Ginsburg functional: 



J{u) 



\Vu\' + G{u) 



dx, where G{z) 



g{s)ds. (52) 







In addition, the Landau-Ginsburg functional satisfies 

dtJ{S{t)uo) = -\\VK{S{t)uo)\\l for < t < T{uo), (53) 
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where K{u) = —uAu + g{u). Furthermore, there exist positive constants /o 
and Co such that — /q < G{s) < Cqs'^ + /o and 



■\\VS{t)uo\\l-fo\n\ < J{uo). 



(54) 



Since J(mo) is bounded below, see (54), and decreasing along orbits, see (53), 
it follows that J{S(t)uo) is a Lyapunov function, see LaSalle and Lefschetz 



(1961). This implies that the mild solution S{t)uo is defined for all t > 0, 
i.e., T{uo) = oo, for every uq G V'^. It is fact that, whenever F = and 
Mo = 0, then uj{uq), the omega limit set of the solution S{t) uq is a nonempty, 
compact, connected invariant set in Qq. 

Let us now return to the CHE with forcing (44). By integrating the 
equation (44), where Uq G Hq and F G 
all t > 0. 



as well, then S{t) uq G 



Hi for 



4.3.2 The Global Attractor 2tn 



Assume that B satisfies the boundary conditions BC, i.e., either the non-flux 



condition (47) holds, or the periodic boundary conditions hold. Let S{t) be 
the semifiow generated by the solution operator on Vq. Then the following 
hold. 

1. S(t) has a nonempty, compact global attractor 2lo in Vq, and 2to at- 
tracts all bounded sets in Vq. The attractor 2lo depends continuously 
on the forcing function F G H^. 



2. When F G 



then the attractor 2to is a compact, invariant set in 



V^"^, for each r with < r < 1. 

The set Qq is nonempty, compact, and invariant with Qo C SIq. 
Lastly, there is an Inertial Manifold for the solutions of the infinite 



(44 


), see 


Foias et al. 


(1988 


) and 


Sell and You 


(2002 



One finds this manifold by using the orthogonal projection P ^ onto the 
lowest N nodes, that is to say, into Span{efc : k < N}, where {cfc} are 
the eigenfunctions for the Bi-harmonic operator B. One then makes 
a change of variables u = v + w, where v = F ^ u, w = Qn u, and 
Qn = I — '^ n- The (f , w) system: 



dtv + Bv = Fn iA{g{v + w)) + F) 
dtw + Bw = Qn {A{g{v + w)) + F) 



(55) 



is equivalent to the CHE (44). One then shows that, for N large, the 



variable w is enslaved to the variable v in some neighborhood of SIq. 
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That is to say, w = for a suitable function $. It turns out the 



the fongtime dynamics of (44) is equivalent to the longtime dynamics 



of the finite-dimensional ordinary differential equation: 



dtV + Bv = Fn iA{g{v + ^(v))) + F) . 



(56) 



See Sell and You (2002) and the references contained therein, for more 
details. 



4.4 BV and EBV in the Higher Nonlinear Regime 

For sufficiently small perturbations and short renormalizing time intervals 
the finite-time Lyapunov vector, the BV algorithm, and the EBV are similar 
in outcomes. As we depart from these conditions, we see significant differ- 
ences in the outcomes of the three methods. The BV outcomes are most 
sensitive to the amplitude and the frequency of the perturbations. Regard- 
less of the amplitude of the perturbations, the EBV outcomes are structurally 
unambiguous and robust. 

We revisit the CY92 simulations of BV and EBV, using the same forc- 



ing and initial conditions used Section |4.2| but set the time scale parameter 
a = 0.01, we increase the spatial resolution to 320 points, and set the in- 
tegration time step at 6t = 0.000015. The perturbation field is the same 



as in (45); however, we will increase the size of the perturbation of each of 



the components. For ej = 0.6,0.8, and 1.2, with j = 1,2, ...,6, in (45) we 
obtain Figure [9} which shows all BVs at t = 19.025. We expect a simple 
structure in the perturbation field and thus the BV should reflect this. The 



EBV results appear in Figure [TOj Comparison of Figures [9] and 10 show how 
the BV outcomes look qualitatively different from their EBV counterparts as 
we increase the initial amplitude of the ensemble. The comparison is further 
aided by reference to Figure [Tl| in which each ensemble member of the EBV, 
at t = 19.025, has been rescaled to 1 in amplitude. The BV outcomes shown 
in Figure |9] do not yield the structural clarity that the EBV ensemble displays 
in Figure [lO] this is a natural consequence of the size-ordering inherent in 
the EBV algorithm. 



In reference to Figure ^ and Figure 10:; we see the similarity between 



one and only one of the BV vectors and the largest EBV. A striking struc- 



tural feature in the EBV results shown in Figure 10 is that there are only 
3 significantly-different shapes, among all of the six vectors (see Figure 11 
where the vector have been rescaled to have size 1 in L^.) For relative com- 
parison, see Figure [9] 

The sensitivity of the outcomes to the size of initial perturbations is 
significantly different in the EBV and the BV outcomes. There are shape 
variations among the BVs as the amplitude of the perturbation is increased. 
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Figure 9: BVs at t — 19.025, corresponding to Cj — e, for all j. (a) e — 0.6; 
(b) e — 0.8, (c) e — 1.2. The outcomes are very sensitive to nonlinearity. 
There is a resulting ambiguity in structure of the perturbation field. 
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Figure 11: EBV outcomes shown in Figure 10. with each vector rescaled to 
an L'^-norm of 1 to aid in visual comparison. Perturbation amplitudes (a) 
e = 0.6, (b) e = 0.8, and (c) e = 1.2. Note that as the perturbation amplitude 
increases, the resemblance between the BV and the EBV outcomes is lost, 
except for one of the vectors. 
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as evidenced in Figure |9] This degree of sensitivity is not as prevalent in the 
EBV outcomes, shown in Figure 11 The BV algorithm does not distinguish 
between different perturbation scales. A large perturbation BV calculation 
will thus not yield information concerning smaller scales. In contrast, see 
Figure 11:;, for the EBV case, where small scale information is evident. (See 
also Figure llOfc). 



5 Implementation Issues 

Like the BV algorithm, the EBV algorithm is capable of dealing with legacy 
code. The important difference between implementing a code that does an 
ensemble of BV and the EBV is that while the former can be run concurrently, 
in the EBV the ensemble members require normalization to each other. In 
terms of coding, this is a minor issue, if the state variable dimension of 
the underlying model is moderate. For very large problems communication 
becomes an issue, but not at all unfamiliar in concurrent or hybrid computing. 
Any additional computational issues borne by the EBV are well outweighed 
by the higher informational content of the EBV over the BV. Moreover, the 
EBV is much more robust under the nonlinear effects, as illustrated in Section 
|4| than BV. This last aspect is very important and it should be studied in 
the context of more chaotic equations in fluid mechanics, meteorology, and 
geophysics. 

The numerical outcome of the BV and EBV algorithms depends on the 



choice of norm used for rescaling. This is alluded-to in Riviere et al. (2008), 
but the reason for this dependence turns out to be easily explained and can 
be significant if nonlinear effects are not negligible. The dependence of the 
outcomes on the norm lies in the fact that it is not possible, in the general 
case, to scale out the norm in the algorithm if the underlying dynamics 
are nonlinear. To illustrate the norm dependence, we consider a simple 2- 
dimensional system 

~df - 
dX 

= -sin(8Xi), t>0, (57) 
subject to the (same) initial conditions (Xi(0), X2(0)) = (0.8,-1). Figure 



12 shows the outcomes of the BV algorithm for three different choices of 
norms, starting from perturbations of size no greater than 0.15. The norms 
we employ are all standard (finite) Lp norms. The resulting BVs reflect 
the characters of these norms. That different outcomes are obtained is not 
surprising: the shape of the unit sphere changes depending on which Lp norm 
is used. 
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(c) 



Figure 12: (a) L2-norm, (h) Li-norm, and (c) Loo-norm of a sample BV as a 



function of time, for the system in (51). The initial conditions are the same 
in all cases, the amplitude of the perturbations was 0.15. 
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A more weighty consideration related to norms concerns physical and 
theoretical considerations. Conservation laws provide guidance for the most 
appropriate norms for the given dynamics (for instance, in Ray leigh- Bernard 
convection, the temperature enjoys a maximum principle, while the velocity 
does not), but other considerations will play a role {e.g., a sup-norm may 
be an obvious choice in determining the location of severe weather events). 
It will not be uncommon, thus, that a mixture of norms may be necessary 
in multi-physics problems; the fact that the choice of norm in the EBV/BV 
affects the results is in fact a good thing, not a bad one. 

To end this section we wish to highlight a practical consideration that 
may be not be familiar to practitioners, implementing BV or EBV algorithms 
computationally. The specific issue is the impact of finite precision computing 
on the outcomes. It is easy to show that these algorithms have high numerical 
sensitivity. In order to illustrate how this plays out we will consider the 
problem of calculating the vectors associated with 

^ = AX, t>0, X{0)=Xo. (58) 

Here A is a square matrix of constants. For linear problems BV, EBV (and 
the finite-time Lyapunov) algorithms must yield the same vectors. We will 
choose to illustrate computationally numerical ill-conditioning on a problem 
that exhibits transient growth due to the non-normal structure of the matrix 
A. We emphasize, however, that the numerical sensitivity we will be high- 
lighting in this example does not hinge on the nature of the dynamics, but 
rather, on the algorithmic form of the BV and EBV themselves. 

We take A to be an upper-triangular Jordan-block matrix of dimension 5, 
where A has a single eigenvalue A, which is repeated on the main diagonal, 
while the diagonal directly above the main diagonal has non-zero entries. 
For simplicity we assume these to have the same value p. We assume further 
that A = — 1 and p = 1. The general solution operator contains "generalized" 
eigen-solutions with growth rates t'^ e~*, where < k < 4. (In this case, the 
eigenvalue has (algebraic) multiplicity 5.) For this linear problem we expect 
to see, provided we take t sufficiently long to forget the transient, a very 
trivial outcome to BV, or EBV. 



In Figure 13 we summarize the results of the BV calculation on this sys- 
tem. We employed an integration time step of 0.001, a perturbation initially 
of magnitude 0.01. We performed the calculation in double precision, using 
an explicit Runge-Kutta 4, but payed little attention to how numerical sen- 



sitivity was handled. In Figure 13 1 we show the BV, i.e., y, as a function of 



time. The calculation of BV, using ([3), clearly diverged, shortly after about 
t = 30 (and thus not shown in Figure 13 1). However, there were indications 



that something was not right even before it became obvious that the solution 
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was wrong: In Figure 13 d we show the 2-norm of M(y„+53^„, 6t) — M{Yn, St). 
It monotonically increases, even though all the factors e~*, for < A; < 4 
in the solution operator, decay for t > 4. 

The outcome is not related to choosing a rescaling time to be too long: 
in fact, in this computation the rescaling is performed at each computational 
step. Moreover, the time step chosen was sufficiently small to guarantee 
asymptotic stability in the numerical integrator. A smaller time step would 
have been able to ameliorate to a certain extent the numerical sensitivity 
of the difference M(F„ + 6yn,St) — M(Yn,St), however, in very large scale 



problems this might not be practical or even possible. Figure 13: shows the 



2-norm of MiYn + 53^n, St) — M{Yn, St), using a numerically well-conditioned 
implementation of the BV algorithm. The strategy used to obtain a well- 
conditioned outcome was to normalize the elements in the required subtrac- 
tion before computing their difference. The well-conditioned calculation was 
capable of qualitatively good results. However, the strategy adapted here to 
increase the numerical stability was by no means generally applicable to all 
problems, nor was it optimal. 



6 Concluding Remarks 

The main thrust of our work is to propose an ensemble-based vector breeding 
algorithm, the Ensemble Bred Vector (EBV) algorithm. It is based on the 



Bred Vector (BV) algorithm introduced by Toth and Kalnay (1993, 1997). 
We compare the EBV to the BV algorithm and the finite-time Lyapunov 
Vector algorithms. In the EBV, an ensemble of initial perturbations is bred 
concurrently and then rescaled by the size of the largest member of the bred 
ensemble. The uniform normalization of all the ensemble members after 
each cycle is the distinctive trait of the EBV algorithm, which leads to some 
profound differences when compared to the BV algorithm. 

As expected, when initial perturbations are sufficiently small, the EBV, 
the BV, and the finite-time Lyapunov vector algorithms lead to similar re- 
sults. We gave a theoretical justification of this phenomenon by looking 
at the corresponding time-continuum analogues of the BV and EBV algo- 
rithms. We rescale frequently: The algorithms and results are formulated, 
for simplicity, assuming that the rescaling is done after every discrete time 
step. 

In Section [3| we develop a solid mathematical basis for both the BV and 
EBV algorithms and show that each algorithm results in good approxima- 
tions of the solutions of the tangent linear model, when the step-size St is 
small. As is seen in Table [T| the EBV algorithm has a substantial advantage 
over the (classical) BV algorithm, in the sense that the drop in the minimal 
EBV error is substantially better than the corresponding drop in the minimal 
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Figure 13: (a) A sample BV, as a function of time. The calculation diverges 
and eventually fails (not shown); (b) 2-norm of the numerically- approximated 
M{Yn + ^y-n) St) — M{Yn, St) as a function of time, for the linear system 
^ = AX. Here, A is an upper-triangular Jordan-block matrix of dimension 
5, where A has a single eigenvalue —1, which is repeated on the main diag- 
onal, while the superdiagonal has entries 1. (c) 2-norm of the numerically- 
approximated M{Yn + 5yn, St) — M{Y„, 5t), well-conditioned case. 
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BV error. 

In the study of the Lorenz attractor, the classical BV algorithm has a 
shortcoming which limits one's ability to use this algorithm to study the dy- 
namics inside the attractor. In particular, equation ([s]) implies that ||53^„+i(t)|| 
e, for all n > and all t G I. The BV algorithm maps the cloud onto a 2- 
dimensional sphere. The third dimension is lost, and with it so are the fractal 
patterns seen in Figures |4] and [Tj Also the many patterns seen in Figure [2] 
are lost because the BV alternative would map everything onto the single 
line at height e. 

Lastly, in Figures [9, 10, and 11, we examine a series of related test prob- 
lems that depart from the tangent linear model. The point here is to get a 
comparison of the performance of the two algorithms, as one moves further 
into nonlinear regime. Once again, one sees a significant advantage of the 
EBV over the BV. 

The theoretical aspects for the EBV and BV algorithms are presented in 
Subsection 2.5. As is noted there, our application to the basic issue of the 
sensitivity with respect to errors in the initial conditions relies heavily on the 
Johnson, et al manuscript, JPS87. 

In conclusion, for the applications described in this article, we believe 
that the new EBV algorithm has been shown to be superior to the traditional 
BV algorithm. Finally, we ask: Is the EBV algorithm the "last" word on 
modifications of the classical BV? Probably not. However, we do expect that 
the EBV will serve as a good starting point for new theories of bred vectors. 
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